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Dynamics of pattern-loaded fermions in bichromatic optical lattices 

Matthew D. Reichl and Erich J. Mueller 
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(Dated: March 24, 2016) 

Motivated by experiments in Munich (M. Schreiber et. al. Science 349, 842), we study the dy¬ 
namics of interacting fermions initially prepared in charge density wave states in one-dimensional 
bichromatic optical lattices. The experiment sees a marked lack of thermalization, which has been 
taken as evidence for an interacting generalization of Anderson localization, dubbed “many-body 
localization”. We model the experiments using an interacting Aubry-Andre model and develop a 
computationally efficient low-density cluster expansion to calculate the even-odd density imbalance 
as a function of interaction strength and potential strength. Our calculations agree with the experi¬ 
mental results and shed light on the phenomena. We also explore a two-dimensional generalization. 
The cluster expansion method we develop should have broad applicability to similar problems in 
non-equilibrium quantum physics. 

PACS numbers: 72.15.Rn, 37.10.Jk, 67.85.-d 


Introduction - An important challenge in many-body 
physics is to understand how interactions and disorder in¬ 
fluence the transport properties of an electron gas. The 
non-interactiM disordered problem was largely solved by 
Anderson [llHI- By studying the expansion dynamics of 
wave packets of weakly interacting atoms, cold atom ex¬ 
periments have found evidence for Anderson localization 
in ID Q and 3D Bi random speckledpotentials and in 
ID quasi random optical superlattices [^ . More recently, 
attention has turned to the interacting problem (?- 25|. 
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FIG. 1: (Color online) Imbalance I = vs time t, 

' ' ■^odd'Y^even 

measured in units of the nearest-neighbor hopping strength J 
for fermions in an incommensurate superlattice of strength A. 
Nodd/even IS the number of fermions on odd/even sites. The 
inset shows the geometry. At time t = 0, I = 1. The dark 
(blue) curves show the result of keeping the first two terms in 
the cluster expansion in Eq. for 20 sites. The light (orange) 
curve shows the result of including three-particle terms in the 
cluster expansion. Red dots correspond to a time-dependent 
DMRG simulation. Here A = 3J, U = 3J, the superlattice 
period = (0.721)“^ and the superlattice phase = 0. 
The density is e = 0.2 in the top graph and e = 0.5 in the 
bottom graph. 


Schreiber et. al devised an ingenious experiment to 
test these ideas. Here we model that experiment. 

The experiment in Ref. [1^ uses lasers to create a one¬ 
dimensional lattice with a weak periodic superlattice that 
is incommensurate with the main lattice (see the inset 
in Fig. [T]). The resulting quasi-periodic potential shares 
features with a disordered one. For example, when the 
potential is sufficiently strong, all single particle states 
are localized. The experimentalists load interacting spin- 
1/2 fermions into some of the odd sites of the lattice, 
leaving the even sites empty. Some odd sites are doubly 
occupied. The atoms hop and interact for time t. The 
experimentalists measure the sublattice imbalance I{t) 

T/,\ Afodd A^even /.. \ 

^ ^ A^odd + iVeven ^ ’ 

where A^odd/even is the number of fermions on odd/even 
sites at time t. In a localized phase, the atoms do not 
travel far from their initial position, and have a rela¬ 
tively high probability of being found at their starting 
point. Consequently in such a phase, one expects I(t) to 
be non-zero at long times. Conversely, in a delocalized 
phase, one might expect I{t) to decay to zero at long 
times. The experiment explores the long time behavior 
of / as a function of superlattice strength and the inter¬ 
action strength. The initial configuration of fermions on 
odd sites is random and the measurements are the result 
of ensemble averages over initial states. The experimen¬ 
talists find two phases: one in which / decays to zero, 
the other in which it is finite. The boundary appears to 
depend on the interactions in a non-monotonic manner. 

In this paper we model the experiment, addressing the 
fundamental question of the interplay of incommensurate 
potentials and interactions. We develop a low-density 
cluster expansion which expresses the ensemble averaged 
imbalance as the sum of terms which involve only single¬ 
particle and two-particle dynamics. Using this compu¬ 
tationally efficient approximation, we numerically calcu¬ 
late the long time imbalance as a function of interaction 
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strength and superlattice strength. Our calculations re¬ 
produce the experimental results and provide insight into 
localization in the interacting system. We also extend our 
method to the case of a two dimensional lattice with an 
incommensurate superlattice in only one direction. The 
extra transverse degrees of freedom give kinetic pathways 
for equilibration; we calculate the consequences. 

Model and Methods - We model the atomic dynam¬ 
ics via the interacting Aubry-Andre model, given by the 
Hamiltonian 14, 2^ 




h.c 


+ A ^ cos(27r/?* -I- + U'^ 


( 2 ) 


The first term describes nearest neighbor tunneling with 
strength J while the second term describes a periodic 
super lattice potential of strength A. For nearly all ir¬ 
rational values of /3, this potential functions as quasi¬ 
random disorder which localizes all single particle states 
for sufficiently large superlattice strength (A/J > 2) 

In this regime, and for infinitely large systems, the sin¬ 
gle particle states are localized with a localization length 
A = (21og^)“^, independent of (3 [13, HI]- If P =plq is 
rational, the eigenstates are extended Bloch waves with 
period q. For large A and large g, the wavefunction in 
each unit cell is sharply peaked, and locally the eigen¬ 
states are similar to the irrational case. 

The localization transition is reflected in the observable 
/ [t ), which for typical irrational /3 and U = 0 relaxes to 0 
for A/ J < 2 but remains finite at long times for A/ J > 2 
(see the inset in Fig. [2]). We define loo = I{t ^ oo). 
Although loo OasA/J 2, the way it vanishes 
depends strongly on /3 and is inconsistent with the naive 
estimate from structureless exponentially localized states 
lest 1/A^ (see Ref. [^, supplementary material). The 
graph of loo vs. P and A/J is fractal (see Fig. SI in the 
Supplementary Information), as it has different behaviors 


J 


for rational and irrational p. Despite this complexity, the 
long time behavior of / is distinct in the localized and de¬ 
localized phase: /(f) captures the localization transition, 
but also probes features of the single-particle wave func¬ 
tions beyond the localization length. 

The third term in Eq. (ED describes on-site interactions 
of strength U. Here we develop a low-density expansion 
to calculate the imbalance in the presence of interactions. 

We define {I[t)) to be the expectation value of the 
imbalance, averaged over the ensemble of initial states, 

1 1 

m) = X-(winz(t)iw) (3) 

{n} 


Here {n} = {fiUi, * 20 ’ 2 , ■■■PnO'n} labels an n-particle ini¬ 
tial state with particles at sites i with spin tr, 
denotes a sum over the ij’s and cj’s, W{{n}) is the 
weight of a given n particle state, Z = 

and hi{t) = e*^‘(iVodd - iVeven)e“*-^* where Nodd/even 
are the number operators (for both spins) on odd/even 
sites. 

To model the experiment, we take W{{n}) = 0 if any 
of the particles are on even sites. We take the initial 
occupation of each odd site to be an independent random 
variable, and hence W{{n}) = e”(l — e)'^'*”", where Ng 
is the number of sites. Our method is readily generalized 
to more sophisticated weights. For instance, as shown in 
Eq. (S12), we can weight the initial states with separate 
probabilities for sites with two atoms (doublons) or one 
atom (singlons) (see also Eig. [3]). 

With this choice of W, the normalization is Z = 1 —(1— 
e)^« which approaches 1 in the iV^ —>■ oo limit. In that 
same limit, the mean density (the number of particles per 
site averaged over the ensemble of initial states) is e . 

Substituting our weight function into Eq. © yields an 
expression for the imbalance as a sum of terms involving 
different numbers of particles: 


m) = ^ ki - E 'c'li} w + 

{ 1 } 


Ns-2 


E'^mW 

{ 2 } 


+ ^(l-e) 


Ns-d 


{3} 




Ns 




{N,} 


( 4 ) 


where C|„}(f) = ({n}|n/(t)|{n}), and the primes on the sums mean they only include odd sites. 

We wish to resum this series, taking advantage of the fact that well-separated particles will move independently. 
Somewhat analogous to cumulants, we define functions C’jn} {t) via 

C{n}it) = Cln}(t) + E ^{1}W+ E C'{ 2 }(t)-|- E C'{3}(i) + •■• + E (- 5 ) 

<{l}eW> <{2}e{r^}) <{3}e{n}> ({„_!}£{„}) 


where X]<{fe}e{n}> denotes a sum over all (^) combi¬ 


nations of k site and spin labels in {n}. We set 
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C{i}(i) = C{i}{t)- These new functions C^k}{t) extract 
the fc-body dynamics from the original functions C^k}(t)- 
First instance, the two particle term C'{iicriy 20 - 2 }(^) = 
cTi, 22 (^ 2 }(^) ^{ 2 i(Ti }(0 ^{ 22 ^ 2 ' 2}(0 thc differencc 

between a term representing the exact dynamics of two 
particles with initial positions and spins iiai and * 2 (T 2 
and the single particle dynamics of a particle initialized 
at site *1 and another particle initialized at site * 2 . In the 
non-interacting limit U = 0, we only have single particle 
dynamics and = 0 for all k > 1. In a diagram¬ 

matic formulation, C involves only connected diagrams. 

Substituting Ep. {|5j) into Ep. Q, and using the argu¬ 
ments in the Supplementary Information gives 

= + ^{2}(^)+0(e^) (6) 

* { 1 } { 2 } 


in the Ng oo limit. For our numerical calculations we 
include the finite size corrections in Ep. (S7). 

Epuation (l6|) expresses the n-particle time dependent 
observable {/(t)) explicitly as the sum of 1-particle terms 
(C'|i}(t)), 2-particle terms (C| 2 }(t)), etc. The first sum 
in Ep. contains Ng terms. The second sum contains 
0{Ng) terms, but when the two particles are farther 
apart than some length scale where ^ is the smaller 
of the one-particle localization length A and the ballistic 
length / = Jt, the particles are effectively non-interacting 
and C'{ 2 } will vanish. Therefore only ^Ng terms con¬ 
tribute to the sum. Similarly, there are only 0{^'^Ng) 
which contribute in the sum over terms. 

Each subsepuent term in Ep. ([6]) is intensive and is 
weighted by a coefficient of the order (the density 
exponentiated to the number of particles in the cluster 
minus 1). This cluster expansion is a non-epuilibrium 
analogue to the virial expansion in statistical physics . 
When the localization length is greater than the system 
size (A > Ng) the series is only guaranteed to converge 
for short times I = Jt < 1/e. Therefore, for calculations 
of the long-time behavior of the imbalance, we focus our 
attention on the localized regime A/J > 2. 

For most of the results in this paper we only keep the 
first two terms in Ep. ([6|). Remarkably, this approxi¬ 
mation, which only involves calculating the dynamics of 
one or two particles, shows all the features seen in the 
experiments of Ref. 261. 

Numerical Results - Figure [1] shows typical {I{t)) for 
interacting fermions in the localized regime. The solid 
blue curves show calculations using the first two terms in 
the cluster expansion in Ep. ®. The imbalance initially 
has a value I{t = 0) = 1, reflecting the fact that the ini¬ 
tial states have particles localized only on odd sites. At 
long times, the imbalance saturates to a non-zero value 
with small fluctuations about the mean. For compari¬ 
son, the red dots show calculations using time-dependent 
density matrix renormalization group (t-DMRG) [s^ ^ |. 
For the DMRG calculations, we average over 100 initial 


states drawn from the probability distribution W{{n}). 
The cluster expansion and the t-DMRG show excellent 
agreement at the smaller density e = 0.2. At the larger 
density e = 0.5 there is some small puantitative disagree¬ 
ment, but the average long-time imbalance is nearly iden¬ 
tical for the two approaches. As a test of the conver¬ 
gence of the cluster expansion, we have also computed 
the contribution from three-particle terms (orange curve 
in Fig. [T]). Including these terms gives small corrections 
to the two-particle calculation and yields better agree¬ 
ment with t-DMRG. 

Figured] shows the long time imbalance /oo as a func¬ 
tion of interaction strength for a series of superlattice 
strengths. We compute (/(t))init by numerically evaluat¬ 
ing the first two terms of Ep. ® at a density e = 0.2. 
Each data point in Fig. d] represents averaged over 

the times 200 < tJ < 500 and averaged over twelve values 
of the superlattice potential phase 4> evenly spaced in the 
range [0,7r]. All simulations were performed on a lattice 
with 20 sites using open boundary conditions. We have 
explicitly verified that finite size effects are negligible; the 
system size was chosen for numerical convenience. 

Each curve is symmetric under U —>■ —U. As pointed 
out in Ref. 32| this symmetry is expected for time- 
reversal invariant operators such as I{t), as long as the 
initial states are localized in space. For \U/J\ < 2A, 
interactions cause some 2-particle scattering states to 
become less localized than 1-particle states, and the 
long time imbalance decreases with increasing interaction 
strength. For larger interactions, the imbalance begins to 
increase again and produces a “W” shape consistent with 
the re-entrant behavior predicted for similar systems . 
The “W” is most pronounced for A/J Ri 3. 

At large interaction strengths, up-spin and down-spin 
particles initially localized at the same site (doublons) 
become bound and have a reduced effective tunneling 
rate JeS ^ J'^/U SQ. The contribution to I^c from 
these doublons causes the long time imbalance at large 
interaction strengths to become greater than the long¬ 
time imbalance at [/ = 0. 

We further explore the contribution of doublons to 
loo by giving doublons and singlons separate weights in 
our average over initial states (see Ep. (S12)). We let 
e be the total density of particles and -q the density of 
doublons. Fig. |3| shows loo as a function of U/J at 
A/J = 3 for three different values of q/e in the initial 
states of the system; 0 (e = 0.5), 0.23 (e = 0.57), and 
0.5 (e = 0.67) for the bottom (blue), middle (orange), 
top (green) graphs, respectively. All other parameters 
are the same as in Figure [H In the case where there are 
no doublons Ioo{U/J = 0) = Ioo{U/J —>■ oo). This is 
a reflection of the fact that the dynamics of singlons in 
the hard core U/ J ^ oo limit is identical to the dynam¬ 
ics of free spinless fermions [l^. As more doublons are 
added to the system, loo at large U/J increases, as ex¬ 
pected from the reduced tunneling rate of bound pairs. 
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FIG. 2: (Color online) Long time density imbalance 7oo as 
a function of interaction strength Uj J for a one-dimensional 
lattice with 20 sites at density e = 0.2. The superlattice 
period is = (0.721)“^ in units of the lattice spacing. The 
different curves correspond to different superlattice strengths: 
A/J = 2, 3,4, 5 (from bottom to top). The inset shows I as 
a function of superlattice strength for U/J = 0. 


The blue and orange points in F ig.pi show corresponding 
experimental results from Ref. |26| . where the doublon 
density was controlled by varying the loading protocol. 

We chose rj and e to best match the experimental data, 
finding excellent agreement. Our best-fit value of rj is 
somewhat smaller than estimates in Ref. [l^. Similar 
discrepancies were seen in DMRG calculations [j^ . 

Motivated by more recent experiments and as a 
further demonstration of our cluster method approach, 
we have extended our calculations to two-dimensional 
lattices. We consider a two-dimensional Hamiltonian 
with a one-dimensional superlattice potential V = 
A cos{2Tr/3ij;+(j)). As before, we take J to be the hopping 
in the x-direction and Jy the hopping in the y-direction. 
In this case we average over initial states where atoms 
are localized on odd sites in the x-direction and are in 
ky = 0 momentum eigenstates in the y-direction. This 
choice of initial states, which requires periodic bound¬ 
ary conditions in the y-direction, was chosen purely for 
numerical simplicity; we expect no qualitative changes 
if we initialize with spatially localized states and use 
open boundary conditions in the y-direction. We once 
again use Eq. ([6]) including only one-particle and two- 
particle terms to compute the even-odd imbalance in the 
x-direction. 

Because the eigenstates are inherently delocalized in 
this situation, we only expect our cluster expansion to be 
accurate for short times. Fig. 2] shows the imbalance I in 
the x-direction, averaged over times between t = 5/J and 
t = 10/J as a function of U/J. These simulations were 
performed on a lattice with 10x10 sites. Scattering in 
the y-direction (transverse to the superlattice potential) 
allows for the density imbalance to relax to smaller val¬ 
ues, and / becomes suppressed as Jy is increased. Similar 


FIG. 3: (Color online) Long time density imbalance 7oo as 
a function of interaction strength U/J for a one-dimensional 
lattice at superlattice strength A/J = 3. The different curves 
show calculations using a cluster expansion on a 20 site lat¬ 
tice with different densities ry of doublons in the ensemble of 
initial states: The bottom (blue), middle (orange), and top 
(green) curves correspond to a ratio of doublons to particles 
of » 7 /e = 0,0.23, 0.5, respectively. The blue and orange points 
are experimental measurements for a small doublon fraction 
(? 7 /e « 0.08) and larger doublon fraction ( 77 /e « 0.5), from 
Fig. 6 of Ref. [^, courtesy of Ulrich Schneider. 



FIG. 4: (Golor online) Density imbalance 7 averaged over 
time from t = 5/J tot = 10/J as a function of interaction 
strength U/J for a two-dimensional lattice with 10x10 sites 
at superlattice strength A/J = 3.0 and density e = 0.2. The 
superlattice potential is only one-dimensional: V{ix,iy) ~ 
Acos(27r/3ia; -|- (/). Jy/J = 0,0.1,1 for the top, middle, and 
bottom (blue, orange, green) curves, respectively. The inset 
shows a diagram of the setup. 


results are observed in Ref. [34| . 

Conclusion - In this paper we have applied a new 
cluster expansion method to simulate experiments 261 


which studied the non-equilibrium dynamics of fermions 
pattern-loaded in quasi-disordered one-dimensional lat¬ 
tices. Our calculations, which involve keeping the first 
two terms in the cluster expansion and account for only 
single particle and two particle dynamics, reproduce all 
experimental features of the long-time density imbalance 
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between even and odd sites, and agree quantitatively 
with simulations using t-DMRG. We have also extended 
our calculations to two-dimensional lattices, finding that 
the density imbalance is suppressed when adding hopping 
in the direction transverse to the superlattice potential. 

Although principally designed to calculate the experi¬ 
mental observable, this cluster approach also gives some 
insight into many-body localization. For example we 
have shown that time dynamics of the many-body wave 
function in the localized phase can be written as a sum 
of 1-body, 2-body, ..., n-body terms. In the dilute limit, 
the dynamics are dominated by few-particle physics, a 
feature which was not previously recognized. 

Our cluster approach can be also used to explicitly 
construct the local integrals of motion which underly the 
^enomenology of the many-body localized phase |15l.ll7l. 
[ 33 , l36j. As detailed below, we use the solution to the j- 
body problem to construct fermionic creation operators 
where = SmTiSra- Our Operators have 

the property that in the f-particle subspace, all of the 

t('j') t(i) tf?) 

OnCT are equivalent for j > i: ana Pi = ana Pi where 

Pi projects into the i particle subspace. Our conserved 

quantities are manifest in the requirement 

[atWaW,F.i7P.] = 0 (7) 


If the al!;'a are “local”, we thereby complete the construc¬ 
tion. 

We take aha^ to create the single-particle eigenstate 
with spin cr and energy suppressing the spin indices 
\n) = |vac). This operator is local if these eigenstates 
are localized. Trivially, Eq. © is satisfied. 

Next we construct 


at(2) = atW + 


( 8 ) 


so that a'lj'a'^Pi = We can always choose the T’s 

such that \na, mr) = aiia'amr |vac) is an eigenstate of H 
with energy Neglecting the spin indices 


^lM = {{j\®{mnl)-6,nhi 


(9) 


There are as many ways of doing this are there are 
ways of assigning the indices to the 2-particle states. We 
choose the indices to maximize the overlap {{n\‘^{m\)\nl). 
If the two-particle states and one-particle states are lo- 
calized, then alia will be localized. Eq. ([7]) is clearly sat¬ 
isfied. Constructing the higher order operators follows 
the same procedure. 


To connect with the existing literature [1^, [ITl, [35 


we note that this construction yields a Hamiltonian of 
the form 


H — ^ ^ ^n'O^na T ^ ^ Unm'ana^nia' T (1^) 

na 'nrn 

crcr' 


where h„cr = limj^oo al^Pana- The coefficients are local, 

meaning -- exp (-max|i„ - i/ 3 |/^k)- They can 

be expressed in terms of the eigenvalues of the fc-body 

( 2 ) ' 

problem; for example Unm = E^^ — Cn — The Sup- 

(T(t' 

plementary Information shows a graph of this quantity 
for typical parameters, illustrating the exponential decay. 
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SUPPLEMENTARY INEORMATION 


Imbalance vs. Superlattice Period in the Non-interacting Limit 


In the non-interacting limit, the experiment is well modeled by the Aubrey-Andre model 

H = - j'^ (4,0-G+i.cr + h.c^ + ^ X] cos(27r/3i -|- 4))c\ ,^Ci^„ (SI) 

i,(7 

where J is the nearest neighbor hopping strength, A is the strength of the periodic superlattice, and j3~^ is the period 
of the snperlattice. As discussed in the main text, this is an interesting model as its behavior depends on if /3 is 
rational or irrational (or in a finite system of length if NgP is an integer or not). 

Starting with a particle on an odd site, we numerical evolve the single-particle wave-function and calculate the 
average long-time imbalance loo = ?^odd — ?T-even, where u-odd/even is the average long-time density on odd and even 
sites, respectively. 

Fig. [ST] shows loo as a function of /3 where Ng = 200, A/J = 3 and (/) = 0. The behavior of the imbalance depends 
strongly on whether /3 is irrational or rational, and thus displays a fractal structure. When Ngj3 = Ngp/q is an integer, 
loo has peaks for even q and troughs for odd q. Increasing Ng leads to finer structure. 
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Derivation of Cluster Expansion 


Here we will derive Eq. (6) given in the main text. From Eq. (4) we have 

2 


(/(t)) e(l - C'{i}(^) +-;r(l - ^ X] ^{2}(^) 


{ 1 } 


{ 2 } 


,3 AT, 

{3} ® {AfJ 

where {n} = {iicri, * 20 ' 2 , labels an n-particle initial state with particles at sites i with spin a, X) {n} 

a sum over the ij's and CTj’s such that the ij’s are restricted to odd sites. 

Substituting Eq. (5) in Eq. (4) we have 


(S2) 


denotes 


Z{I{t))—e{l — C{i}(t) + y (1 - |^C{ 2}(0 + X/ 

{1} {2} ({i}e{2}> 

'X! [^{3}W+ X] ^{2}W+ X ^{i}W 

{3} <{2}G{3}) <{1}G{3}) 




(S3) 


where X)({fc}G{n}) denotes a sum over all ()() combinations of k site and spin labels in {n}. For example, neglecting 
spin indices: X;{ 2 }' E{{i}e{ 2 }> /({!}) = Enodd(/(H) + fih))- 

i2odd 

We note the following identity: 



(S4) 


where the combinatorial factor is the number of ways of choosing the n — k elements of {n} which are not in {fc} out 
of the Ns — k available starting positions/spins. 

Substituting this identity into Eq. (S3) yields 


Z{I{t)) =e(l-e) 


£3 


^X ^[X ^{2}(i) + f ) X 

{1} {2} ^ '' {1} 

[X'^131(0 + ") X'^mW + (\') X'^mW 

{3} ^ ^ {2} ^ ^ {1} 


+ ... 


Collecting like terms, we have 


which can be expressed as 


^w))=f: h"(i - d~--”("‘_.') E'cuiu) 

n=l ^ ^ {1} 

n=2 ^ 2 ^2} 

n=3 ^ 2 _j- 3 j 


+ ... 


(S5) 


(S6) 


(/(t)) — Hi(e) X C'{i}(t) + ^2(6) X C'{2}(0 + 213(e) X ^{ 3}(0 + "-+X 2lAfs(e)C'{Arj(<) 

{1} {2} {3} {N,} 


(S7) 








where Ak{e) = Y^n=k Taking the Ng ^ oo limit gives Eq. (6) to O(e^). Including finite size 

corrections, we have Ai{e) = and ^ 2 (e) = ■ 


Doublon Weighting 

Here we develop a cluster expansion for an ensemble averaged imbalance (lit))' which weights initial states with 
separate probabilities for doublons and singlons. We define by 


(Ht))' = 


1 


{N)Z 


E E'l 


{n-f}{n^}{n^}P 


rid 


{{nt}{ni}{nd}\ni{t)\{n^}{rn}{nd}) 


(S8) 


n-^-\-ni-\-nd<Ns/2 


where {n-f} = {ii,i 2 i-An^}, {*^ 4 ,} = {jhh,- JnJ, and {n^} = {fci, fca,label the sites of up spin singlons, 
down spin singlons, and doublons, respectively. The symbol ^ ^ possible locations 

of up-spin singlons, ni down-spin singlons, and Ud doublons, restricted to odd sites, p and r are weights for the 
singlons and doublons. Z is a normalization factor given by 


Z = 


E 


n^,ni,nd 

n-^+ni-\-nd<Ns/2 


pnt+n^^ridf \ ^ 2p-f 

yn^niUdJ 


m 


^liere IS the number of ways of assigning n^ + n^ singlons and Ud doublons to 7V,/2 

odd sites. (N) is the mean number of particles and is given by 


{N) = 


Z 


E 




n-f,rn,nd 
nf+n^,+nd<Ng/2 


f Ns/2 \, , 1 

]{n^ + ni + 2nd) = 7 
yn^niudj 1 


Ns{p + t) 
+ 2p + T 


(SIO) 


We define C!^„i}{ni}{nd}(.^) = {{''T't}i''^i}i''^d}\ni{t)\{nf}{n^}{nd})■ We decompose the expectation value C{t) into 
single particle contributions, two particle contributions, etc. in a manner similar to Eq. (5) in the main text: 


= E C'{i}{o}{o}(t) + E ^{o}{i}{o}(t) 
({1}6W}) ({1}6W}) 

+ E C'{o}{o}{i}(t) + E ^{i}{i}{o}(t) + ■•■ 
<{i}e{"4> ({i}e{"t}di}eW}> 


(Sll) 


denotes a sum over all labels in {k} and C'{i}{i}{o}(t) = C'{i}{i}{o}(t) - C'{i}{o}{o}(t) “ C'{o}{i}{o}(t)- 
There are higher particle number terms in this decomposition, but for the low density limit we consider here, it is 
sufficient (and notationally simpler) to keep terms up to two-body. We note that two-body terms like C'{2}{o}{o}(0 = 
C{ 2 }{Q}{o}(i) ~ ^{i}{o}{o}(0 “ ^{i}{o}{o}(i) vanish, since two atoms with the same spin do not interact. 

Substituting Eq. (ISlip into Eq. (IS8I) and performing simple summations yields 


int))' = 


1 


-f 


p -I- T TV, 

V. -E' ’C{o}{o}{i}(i) 


E '‘^{iHolfo} W + E '‘^Wfiim W 

{1} {1} 

P + T 1 


(S12) 


p -I- r TVs 


{ 1 } 


(1 -f 2p -f r) TVs 


E{i}di}^{iHi}{o}(^) 


We vary p and r in Eq. (IS12P to produce Eig. 3 in the main text. 


Local Integrals of Motion 


( 2 ) 

Eig. S2 shows the coefficients U/nn 

n 


which appear in Eq. (10) of the main text. 
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FIG. S2: Two particle interaction term appearing in Eq. (10) of the main text. Here A/J = 3, U/J = 3, and /3 = 0.721. 

Darker colors correspond to larger values of \Um}i\- For large |n — m|, is exponentially small. For n = m, Uml U. 




